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APPLICATIONS OF NONVARIATIONAL FINITE ELEMENT 
METHODS TO MONGE-AMPERE TYPE EQUATIONS 

TRISTAN PRYER 



Abstract. The goal of this work is to illustrate the application of the nonvari- 
ational finite element method to a specific Monge- Ampere type nonlinear par- 
tial differential equation. The equation we consider is that of prescribed Gauss 
curvature although the method can be generalised to any Monge-Ampere op- 
erator. 

1. Introduction and problem setting 

The nonvariational finite element method (NVFEM) introduced in |LPllaj is a 
numerical method aimed at problems of the form 



Aixy.BMx) = fix) (1.1) 

where for each x e f2 the matrix A{x) G Sym([R^^^), the space of bounded sym- 
metric positive definite matrices and where D'^u denotes the Hessian of the function 
u. The operation B.C = trace {B^ C) is the Frobenius inner product between two 
J2 d X d matrices. Classical finite element methods are applicable to this problem if 

^ we assume the coefficient matrix A is different iable. In this case we may rewrite 



^^ (1.1) in variational or divergence form via the introduction of an advection term 

^_^ since 

> f{x) = A{x):D^u{x) = div {{A{x)\/u{x))) - DA{x)\/u{x) (1.2) 

^J^ where 

^ / d d \ 

^ BA{x) = I ^diai^i{x), . . . ,^diai^d{x) I . (1.3) 

^^ Note that we are using the convention that \/(j) = (dicj)^ . . . ^dd(j)y is the column 

^Sj vector formed of first order partial derivatives of a ^-multivariate function (j). 

i-H The introduction of the advection term may result in the variational problem be- 

^ coming advection dominated. This is undesirable in the finite element context and 

• ^ stabilisation terms become necessary to derive a viable numerical method [EG041 

KN c.f.]. Interestingly if ||DA||l_(^) > ||A||l^ ^ applying the NVFEM to Ol) does 

C^ not result in an unstable scheme, whereas applying a standard FEM to ( |1.2p does. 

This is numerically demonstrated in [LPlla, §4.2]. It may even be the case that A 
is not different iable, in which case the standard FEM cannot be applied. 
The fully nonlinear problem 

^(D^ii) = (1.4) 

is related to the nonvariational problem ( |1.1[ ) by the fundamental theorem of cal- 
culus. If ^ is sufficiently regular and u solves (|1.4|) then u also solves 



/ ^'{tTi'^u)dt 
Jo 



:D^ii + ^(0) = 0, (1.5) 
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where ^' is the Frechet derivative of ^. It is also the case that a Newton hneari- 
sation apphed to ( |1.4| ) results in a sequence of linear nonvariational PDEs: 

^'(D2ii^):D2 (2i^+^ - u"") = -^(D^ix^). (1.6) 

The Monge- Ampere operators are an extremely interesting class of fully nonlinear 
PDE. These arise from differential geometry and optimal transport problems; they 
take the form 

^(D^u) := det {D^u) - f (Vix, u, x) = 0. (1.7) 

For example the Monge-Ampere-Dirichlet (MAD) problem is the case when / = 



f{x) and (1.7) is coupled with a Dirichlet type boundary condition {u = g on dQ). 
This particular equation is a prototypical example of a fully nonlinear PDE. 

There are a variety of numerical methods available for the more general Monge- 
Ampere class of fully nonlinear PDE \1.1\ . In |0beQ 8' the author proposes a wide 
stencil finite difference scheme. In [BohOSj a C finite element scheme based on 
the Argyris element is used. In a series of papers Feng and Neilan |FNQ9b[ IFNQQaj 
construct numerical approximations of solutions to sequences of quasilinear bihar- 
monic equations. This is very reminiscent of the vanishing viscosity method first 
studied for use in fully nonlinear first order PDEs. The method is aptly named the 
vanishing moment method. More recently in [ BGNSll] a consistent penalisation 
method has been introduced for these problems. Finally, Awanou [Awall uses a 
Laplacian relaxation technique to study these equations. 



For the Monge- Ampere type equation (1.7) to be well posed we require i? C R^ 



to be a convex domain and / > 0. The Monge- Ampere operator is elliptic over the 
cone of strictly convex functions in f2 and under the constraints above will admit 
a unique convex viscosity solution [CC95 . 

In this work we will study the equation of prescribed Gauss curvature. This 
arises from the problem of finding a function u such that the graph of u has a 
specified Gaussian curvature K. In this case we have that 

K= —-- 1.8 

and hence 

^{D^u, Vu, x) := det "D^u - K (l ^ \^u\^) • (1-9) 

Note that K = K{u,x). 

The linearisation of this problem can be calculated in a direction v as 

^'(D'^u, Vu, xy.'D'^v = lim - (^(D'^u + eD^v, \/u + eVv, x) - ^(D'^u, Vu, x)) 
= cof D^uiD^w + {d + 2)K({l + |Vw|^] (Vw)'^Vw 

a.io) 

and thus the linearisation is elliptic if coiD^u is an elliptic operator. This holds for 
convex u. 

2. Discretisation 

The process of discretisation can be sought in two ways. We may look at the 
continuous problem and discretise this directly, resulting in a system of nonlinear 
equations, or we may first linearise the problem and discretise from there. Dis- 
cretising the nonlinear problem directly is certainly possible but is more technical, 
as discussed in |BGNSll] . For brevity we will perform a Newton linearisation on 
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1.7|) and discretise the sequence of linear nonvariational PDEs in a similar light to 



ILPllbj . 

Let ^ be a conforming, shape regular triangulation of i7, namely, ^ is a finite 
family of sets such that 

(1) K G 3^ implies K is an open simplex (segment for d = 1, triangle for d = 2, 
tetrahedron for d = 3), 

(2) for any iC, J G ^ we have that i^ H J is a full sub-simplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) [j^^^K = 72. 

We also define (o to be the skeleton of the triangulation, that is the set of sub- 
simplices of ^ contained in f2 but not df2. For d = 2, for example, ^ would 
consist of the set of edges of ^ not on the boundary. 

We use the convention where h : f2 ^ R denotes the meshsize function of 5^, 
i.e., 

h{x) := mdiKhK- (2-1) 

K3x 

2.1. Definition (FE spaces). Let P^(^) denote the space of piecewise polyno- 
mials of degree k over the triangulation 5^ of i7. We introduce the finite element 
spaces 

V = FP{^) n C^{f2) n RliQ) and W = P^(^) n C\f2) (2.2) 

to be the usual space of continuous piecewise polynomial functions and 

S:=VxW'^'. (2.3) 

2.2. Remark (generalised Hessian). Given a function v G H^(i7), let n : df2 -^ 
R^ be the outward pointing normal of i? then the Hessian of 'U, D^'u, satisfies the 
following identity: 

(D^^lc^) = - / Vv(g)V(^+ / Vv(g)n(/) V(^GH^(i7) (2.4) 

Jn JdQ 

where (• | •) denotes an appropriate duality pairing. It follows that we may weaken 
the regularity assumptions to 'U G H^(i?) fl H^(9i7). 



2.3. Definition (finite element Hessian). From Remark 2.2 and in view of Reisz 
representation theorem we may define the finite element Hessian such that 



/ 



H[V]^ = {'D'^V\^) V<I>gW. (2.5) 

2.4. Proposition (symmetry of the finite element Hessian). The finite element 
Hessian is symmetric, that is for each V eV 

[ H[v]^ = [ {H[v]y^ y^ ey. (2.6) 

In view of the constraints to the continuous problem (1.9) to admit a unique 
solution it is also necessary to construct a discrete notion of convexity. This has been 
developed in |AMQ9j and is naturally passed down from the concept of distributional 
convexity. 

2.5. Definition (finite element convexity |AMQ9j ). A function, v G H^(i7) H 
H (9i7), is said to be finite element convex if 



/. 



H[v]^ is positive semidefinite V ^ G W (2.7) 



where $ > on i?. It is strictly finite element convex if (2.7) is positive definite. 
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Definition 2.3 allows us to construct what is essentially a 2-0 mixed method 



where the Hessian of the solution to (1.7) is treated as an auxiliary variable in the 



formulation (as opposed to the 1-1 mixed methods commonly found in the literature 
by decoupling a second order PDE into a system of first order PDEs [BF91, c.f.]). 
Given the linearisation ( |1.1Q[ ) we formulate the problem in the discrete setting as 
follows: Given an initial guess (/7^, HlU^fj G S that is strictly finite element convex 



(2.7), forn G N find {U'',H[U'']) e S such that 



{H[U''],^) + (V/7^ V<^) - Y^ (V/7^ ^n^)^ = (2.8) 

eedf2 

(2.9) 



Where for the problem of prescribed Gaussian curvature ( |2.9[ ) is 

= / cofif[C/"-^]:if[f/"-C/"-^]* (2.10) 

Jn 

+ f 2dK(l + {VU^-^^Y^^ (Vf/"-i)V ([/" - C/"-i) $ 



/ 



/ 9\ (d+2)/2\ 

detff [/7^-^] -Kh^ |V/7^-Y) ) ^• 



Due to the symmetry property given in Proposition |2.4| we may simplify the 
problem somewhat to seeking only the upper (or lower) triangular parts of the 
finite element Hessian. This reduces S = V x W*^^ -^d)/2^ 

2.6. Theorem (solvability of the discrete system |PrylQ| ). Let [/ G V be the 
nonvariational finite element approximation to u, the solution of 

A-.B^u = f, (2.11) 

where A is an elliptic operator. Then we have a discrete inf-sup condition, that is 
the linear system is always invertible. Hence, assuming the linearisation maintains 



ellipticity the discrete problem (2.8)- (2.9) is well posed. 



3. Numerical experiments 



In this section we detail numerical experiments on the formulation (2.8)- (2.9). 

We will consider the case d = 2 and when K > is some prescribed curvature. 
In each of the experiments we choose p = 2, i.e., V consists of piecewise quadratic 
functions. The domain i? is taken as a square whose size differs on each of the 
experiments and the triangulation ^ is unstructured. All of the numerical exper- 
iments have been conducted using the DOLFIN environment of the finite element 
package FEniCS [LWIO . 



In Figures |l||2| we construct classical solutions to (1.9) in order to look at the 
numerical convergence of the method. In Figure [3] we consider i^ as a constant over 
the domain [—.57, .57]^. These results can then be compared with the two other 
numerical studies found in the literature on prescribed Gauss curvature |FNQ9b[ 



lAwallj . In these experiments the authors note that the problem (1.9) is well posed 
only for K < K^^^ and estimate the value of K^^^ by asserting when the numerical 
algorithm proposed breaks down. 

The initial guess to any Newton iteration is paramount due to the well known 
overshoot property. In the case of Monge- Ampere type linearisations it's especially 
important since (discrete) convexity must be maintained during the iterative pro- 
cedure for the problem to remain well posed. In each of the tests below we initialise 
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the algorithm by approximating the solution of the MAD problem over the initial 
mesh as detailed in |LPllb| . 



detD'^u = K in i? 
u = a on df2. 



(3.1) 



The iterative procedure given by the discrete problem (2.8)- (2.9) is terminated 
when two concurrent iterates satisfy ||[/^ — /7"^~"'^|L ,^, < 10""*^^ 



Figure 1. In this experiment we fix choose a convex solution u 
whic h classically solves the equation of prescribed Gauss curvature 
(1.9) over the square [—.5, .5]^. That is, we fix i^ = \x\ and calcu- 
late K = K{x^u). We solve the discrete problem over a sequence 
of concurrently refined meshes and ascertain the errors and conver- 
gence rates for the problem in L2(i7), Ho(i7) and a discrete H (i?) 
seminorm. Notice that ||i^-[/^|| ^ 0{h^), \u-U^\^ ^ 0{K^) 
and p''u-H[U^]\\ ^0{h^-^). 
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(a) Errors and convergence rates for the prob- 
lem, p = 2. 



(b) Solution plot 
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Figure 2. In this experiment we fix choose a convex solution u 
whic h classicahy solves the equation of prescribed Gauss curvature 
(1.9) over the square [—.5, .5]^. That is, we fix li = exp (\x\ /2j 

and calculate K = K{x^u). We solve the discrete problem over 
a sequence of concurrently refined meshes and ascertain the er- 
rors and convergence rates for the problem in L2(i?), Hj(i7) and 



a discrete H (i?) seminorm. Notice that 



U 



N\ 



0(h^) and \\B''u-H[U^]\\ ^0{h^-^). 
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(a) Errors and convergence rates for the prob- 
lem, p = 2. 



(b) Solution plot 
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Figure 3. In this experiment we fix /i ^ 0.009 and imple- 
ment the discrete problem over an unstructured mesh of the square 
[—0.57,0.57]^ consider the case K is constant. We choose various 
values of i^ > and display a contour plot together with a side 
view of the discrete solution. Note that the numerical algorithm 
fails to converge for i^ = 2. 
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(a) Contour plot for K = 0.01 




(b) Solution plot for K = 0.01 
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(c) Contour plot for K = 0.1 



(d) Solution plot for K = 0.1 
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(e) Contour plot for X = 0.5 




(f) Solution plot for K = 0.5 
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Contour plot for X = 1.0 



(h) Solution plot for K = 1.0 




(i) Contour plot for X = 1.5 




(j) Solution plot for K = 1.5 



